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Abstract. In this paper we describe the equations of motion developed for a point-mass 
zero-thrust (gliding) aircraft model operating in an environment of spatially varying atmospheric 
winds. The wind effects are included as an integral part of the flight dynamics equations, and the 
model is controlled through the three aerodynamic control angles. Formulas for the aerodynamic 
coefficients for this model are constructed to include the effects of several different aspects 
contributing to the aerodynamic performance of the vehicle. Characteristic parameter values 
of the model are compared with those found in a different set of small glider simulations. We 
execute a set of example problems which solve the glider dynamics equations to find the aircraft 
trajectory given specified control inputs. The ambient wind conditions and glider characteristics 
are varied to compare the simulation results under these different circumstances. 

Section 1: Introduction 

Small unpiloted aircraft (UAVs) are conceptually attractive for surveillance, communication 
relay, and other loitering-dominated missions. This is because the lack of an on-board pilot 
eliminates risk to human life, and because vehicle cost tends to shrink with the size of the 
vehicle. In the missions considered for such aircraft, a substantial capacity for loitering endurance 
is needed and, indeed, it can be asserted that there is no practical upper limit in the amount of 
loiter capacity that is desirable. 

Small UAVs currently contemplated or deployed typically use propulsion to provide energy 
to the airframe. This requires fuel to be carried (placing a ceiling on endurance and driving 
system mass up) or requires that the vehicle carries a “sufficiently large” array of solar cells and 
batteries to stay aloft. Difficulties with these options motivate exploration of alternate paths to 
mission endurance. 

One of these is to extract energy from the ambient atmospheric velocity field through which 
the aircraft flies, as is done naturally by large soaring birds - consider, for example, vultures 
thermalling over carrion. This option is attractive because, superficially, the energy is “free” 
- assuming that the aircraft is configured to soar well (like a sailplane), it need only fly in a 
certain way and it will gather sufficient energy from ambient atmospheric motion to stay aloft. 

The complication encountered here is that the aircraft’s flight path must simultaneously 
satisfy mission performance requirements and atmospheric energy extraction requirements. A 
simulation model of such a vehicle is needed in order to carry out the analyses leading to such 
mission designs, and one is provided in this paper, for the case of purely gliding flight. The key 
features of the model provided are as follows: 

• The airframe dynamics are represented in point mass form, fully accounting for the three- 
dimensional velocity of the atmosphere, an arbitrary function of location in an Earth-fixed 
coordinate system. The atmosphere representation can be time-dependent, if desired, but 
since the equations of motion are autonomous, such time dependence will not appear 
explicitly. 

• A parametric model for aerodynamic lift, drag, and side force coefficients is constructed, 
based on a similar model of very small aircraft [1]. The model is parameterized on vehicle 
wingspan (between 60 and 140 inches) and aspect ratio (between 6 and 16). The weight 
of the vehicle is also used as a variable parameter in the overall glider model. 
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These features differentiate the above model from those typically provided for powered flight of 
full-sized vehicles. In the latter case, the vehicle’s velocity is typically overwhelmingly larger 
than that of the atmospheric wind through which it flies; winds, in this case, can be either 
ignored or treated as small disturbances. For small gliders, on the other hand, the norm of the 
ambient wind may quite possibly exceed that of the vehicle’s velocity. Similarly, it is important 
that the model reflect the aerodynamics of small, slow, vehicles given that drag in particular is 
a strong function of Reynolds number. 

The glider model will be described in the next three sections of the paper, beginning with the 
basic equations of motion in Section 2. The basic aerodynamic force equations (without wind) 
are then described in Section 3, followed by the addition of wind influence and construction 
of the full model in Section 4. In Section 5 the method of finding values for the aerodynamic 
force coefficients is described and some resulting characteristic parameters are compared with 
a literature source. The behavior of the glider model is demonstrated in Section 6 with a set 
of example problem results looking at solutions of the flight dynamics equations with fixed 
control inputs. Section 7 gives some overall conclusions about the work described in this paper. 
Appendix I contains more complete details of the glider model. 

Section 2: Glider Equations of Motion 


The standard form equations of motion for this point-mass zero-thrust aircraft, given a flat 
Earth assumption, are [2]: 


h = 

V sin 7 

X = 

V cos x cos 7 

y = 

V sin x cos 7 

v = 

D 

g sin 7 

m 


7 = — — (L cos o + C sin a) — -f- cosy 

mV V 

y = (L sin a — Ceos a). 

mV cos 7 


( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 

( 6 ) 


These equations define the inertial motion of the aircraft, with g being gravitational acceleration, 
m the vehicle mass (a known glider characteristic), V = ||F|| the vehicle speed, and D , C and L 
the aerodynamic force components (functions of the glider speed and angles as will be discussed 
below). 

Equations (1-3) give the inertial position rates in an Earth-based coordinate system, also 
representing components of the velocity V = [x y z] T in the directions of the axes in the unit 
vector set 
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(north, east, and downward respectively), with the modification that instead of the downward 
£ component, equation (1) represents the change in the altitude h = —z. Equations (4-6) define 
the vehicle dynamics, specifying the rates of change of the glider speed and direction. The 
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azimuth angle x and elevation angle 7 , along with the bank angle a (a control rotation about 
the velocity vector V), form the euler sequence {x, 7 , c}. This moves from the inertial frame 
axes to the right-handed set of “velocity frame” axes Uy = \x v y v iy], which are aligned with 
the inertial velocity vector so that x v = Vy||F||. This rotation is represented by: 

Uy = Ryi(x, 7 , (7 )Uj = Ryj (8) 

(since we are in inertial coordinates and so as noted above U/ = I). Details of this rotation as 
well as the rest of the equations of motion formulation are given in Appendix I. 

Section 3: Aerodynamic Forces With No Wind 

The drag, side, and lift forces D , C , and L make up the separate velocity frame components 
of the overall aerodynamic force vector F \ v = —Dx v — Cy v — Lz v . Control of the vehicle model 
takes place through modulation of the aerodynamic force vector and by rotation of the C and 
L components via the bank angle a in equations (5-6). D , C , and L are functions of the 
aerodynamic angles a (angle of attack) and (3 (sideslip angle), as well as the glider speed V and 
a few aircraft and atmospheric parameters. 

Along with a, the angles a and [3 assume the role of control variables, by changing the 
direction from which the airmass impacts the vehicle body. The sequence of angles {—/3,a, 0} 
is used to rotate (about the velocity axes) from the velocity frame axes to the vehicle “body 
frame” axes Ub = [27 jjb ' z b\- This rotation is represented by: 

Ub = (TivR B v(-/3,a,0)Tfj) U v . (9) 

Note that the j3 rotation is negative; this is a convention which we will follow in order to preserve 
consistency with the flight dynamics literature [2], The coordinate transformation matrix Tjy 
(which changes a vector from velocity coordinates to inertial coordinates) and its inverse in 
equation (9) are necessary since the —j3 and a rotations are performed about the velocity axes 
but we begin and end in inertial coordinates. The transformation matrix is in fact given by 
Tjy = Uy = Ryi , so that equation (9) is rewritten as 

Ub = ( RyiRsvRy] ) Uv- (10) 

In the absence of wind, the drag, side, and lift forces in the velocity frame are given directly 
by the chosen formulas 


D 

= qSC D (a, (3) 

(11) 

C 

= qSCc((3 ) 

(12) 

L 

= qSC L (a ), 

(13) 

where S is the glider wing surface area, 

and the dynamic pressure q is given by 




(14) 


with atmospheric density p. These force values then can be used in the system (1-6) to express 
the glider flight dynamics in zero-wind conditions. 
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The aerodynamic coefficients Cd, Cc, and Cl are functions of the aerodynamic angles in 
a manner to be discussed in Section 5. This formulation includes more detail than previous 
models of glider flight in wind [3, 4, 5] which use the coefficient Cl (and Cd as a function of it) 
directly, rather than as a function of the aerodynamic angles. This extension to angle of attack 
and sideslip angle dependence will create additional complexity in the present model when wind 
becomes involved, as will be seen in the next section. 

Section 4: Aerodynamic Forces With Ambient Wind 

The presence of ambient wind changes our physical interpretation of a and (3. They remain 
control variables with a, but rather than defining the incidence of the body axes to the incoming 
airmass, they now only describe the attitude of the body with respect to V, and will be labeled 
the “velocity-relative” control angles. The definition of the actual attack, side, and bank angles is 
in reference to the incoming airflow direction, and so with the addition of ambient wind altering 
the direction from which the air impacts the vehicle, their meaning will diverge from the velocity- 
frame-based definition in the last section. We introduce a “wind-relative” reference frame to 
serve as the midway point in a transformation from inertial to body coordinates that includes the 
wind effects. Wind-relative force components D Wl C w and L w must be computed using newly- 
found wind-relative angles, before the results are transformed back to velocity coordinates for 
use in the equations of motion. 

Denoting the ambient wind vector as w(h,x,y), the wind- relative glider velocity (velocity 
with respect to the moving airmass) is V w = V — w. We will need to find an euler rotation 
sequence {Xw,lwi &w} which takes the inertial frame axes to the wind-relative frame axes Uw = 
[x w Vw z w \ (with x w equal to 14;/||14 j||), as well as a sequence {—(3 Wl a w . 0 } in the wind-relative 
frame which takes the wind- relative axes to the vehicle body axes. Equation (14) will use the 
wind-relative speed V w = ||t4,|| instead of V, and equations (11-13) will use a w and (3 W to find 
D w , C w and L w . These rotational sequences, parallels to equations ( 8 ) and (10), are written as: 


U\y — RwiiXwiIwi &w)Ui, 

(15) 

Ub = [RwiRbw{— P% iPi>Ot v , 0 )R\yi ) Uw- 

(16) 

The angles Xw and 7 ,,, are found by ensuring that x w = RwiX . The 

other wind-relative 

angles ( a w , a w and (3 W ) must be found by comparing the transitions to the 
velocity frame and the wind-relative frame, so that: 

body frame via the 

(RwiRbwRwi) R-WI = (RviRbvRvi) Rvh 

(17) 

which reduces to 


RwiRbw = RviRbv. 

(18) 

Details of the solution for the angles are given in Appendix I. 


The drag, side, and lift forces are calculated starting with equations (11-13), now in terms 

of the wind-relative angles a w and /!,„ , as: 


Dw = qSC D (a w ,/3 w ) 

(19) 

C w = qSCc(Pw) 

( 20 ) 

Rw — qSC L (a w ). 

( 21 ) 
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( 22 ) 


The dynamic pressure q is given as a function of the wind-relative speed by 

1 Tr2 

(1= 2^- 

However, we want the forces in terms of component values along the velocity axes (F\ v = 
[-D —C — L] T ). Therefore we apply a transformation from F\ w = [—D w —C w — L W ] T , in 
wind-relative coordinates, through inertial coordinates to the final form of the forces: 

Fa v = RyjRwiFA w - (23) 

These velocity frame force components are substituted into equations (4-6), along with the 
original velocity-relative bank angle a, to control the vehicle dynamics. 

Section 5: Aerodynamic Force Coefficient Estimation 

The aerodynamic coefficient model is similar to that of the simulation-based soaring perfor- 
mance study of R/C (radio-controlled) model sailplanes by Beron-Rawdon [1], The model here 
uses values of parameters from Beron-Rawdon’s paper, and the following formulas for the con- 
struction of aerodynamic coefficients include most of the constitutive elements of that model as 
well, although the two are not identical. The representation of aerodynamic effects contributing 
to the following formulas can be found in sources such as Etkin [2], Perkins and Hage [6], and 
Anderson [7]. Airfoil performance data is taken from a low Reynolds number airfoil study by 
Selig et al [8] on which Beron-Rawdon’s models were also based. The drag, lift, and side forces 
will be functions of the angle of attack and side angle of the aircraft (a and /3, or when including 
wind, a w and f3 w ), and certain design parameters. 

The main design parameters are wingspan l (in), aspect ratio AR , and weight W (lb), which 
can be varied to test differently structured gliders. Beron-Rawdon uses ranges of i from 60 to 
140 in, AR from 6 to 20, and W from roughly 0.5 to 12 lb. The following values are also chosen 
for simulated small gliders: 

• fuselage area Sp (in 2 ) (a range of values as a function of wingspan, specifically 86, 145, 
216, 300, and 396 in 2 for wingspan values of 60, 80, 100, 120, and 140 in 2 respectively), 

• fuselage moment arm length (from vehicle center of gravity to both horizontal and vertical 
tail aerodynamic centers) = 0.28£ (in), 

• horizontal tail volume ratio Vp = 0.4, vertical tail volume ratio Vy = 0.02, 

• mean aerodynamic chord c = 1.03 l/AR (in), 

• Oswald efficiency factor e = 0.95, 

• “infinite-span wing” airfoil aerodynamic characteristics from the SD7037 airfoil [8], which 
for Reynolds number 150,000 has values of lift curve slope «o = 0.1(180/7r) (1/rad), zero 
point a 0 = -2.5(7r/180) rad, and wing profile drag given by C d = C do +C dL {C L -C Lmin ) 2 = 
0.01 + 0.05(C l - 0.4) 2 , 

• fuselage drag coefficient Cp> F = 0.008, tail drag coefficient Cp> T = 0.01, and miscellaneous 
“extra drag” coefficient Cp> E = 0.002. 
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Initial calculations find the area values for the various surfaces from the given parameters. 
The wing surface area is S = l 2 /AR (in 2 ), the horizontal tail surface area is St = V d jcS/it (in 2 ), 
and the vertical tail surface area is Sy = VylSjlt (in 2 ). 

Prom the “infinite-span wing” airfoil lift slope value a o, the finite wing slope is 

^ = 1+ aJleAR) (24) 

(in units of 1/rad). The angle of attack value where the finite wing lift coefficient 67 ,( 0 -) = 0 is 
the same value as for the infinite-span wing, «o. Thus the glider lift coefficient is given by the 
following function of or. 

C L (a) = C Lq (a - tt 0 ). (25) 

The drag coefficient Cd{ol) is built up of several elements. In the vertical plane (ignoring 
side angle effects for now), these are as follows: 

• constant- valued fuselage, tail, and “extra” profile drag values given by Cd f S f / S, Cd t (S f + 
Sy) / S , and C Fe (the multiplications by surface area ratios are due to the fuselage and 
tail drag effects acting on those smaller surfaces, not the wing area, and so the use of the 
ratios will effectively replace the wing area in equations (I I- 13) by the fuselage and tail 
areas for the relevant terms), 

• wing profile drag C d = C do + C dL ( C L - C Lmin ) 2 , from the airfoil data, 

• induced drag from lift of G\) n = 6|/ (neAR). 

Added together these form the total drag from constant and a-related influence. The constant 
part of the formula is 


C Do = C Df S f /S + C Dt (St + Sv)/S + C De + C do , (26) 

with the overall formula given by 

C D = C Do + C dL (Cl - C Lmin ) 2 + C 2 L /(ireAR). (27) 

When the side angle j3 is nonzero, there will also be a side force and additional drag, pri- 
marily resulting from the vertical tail acting like a wing, providing “lift” in a sideways direction. 
Side angle effects are not considered in Beron-Rawdon’s paper, so we extrapolate here a roughly 
approximate formulation. The aspect ratio of the vertical tail is usually much smaller than that 
of the wing, and here is assumed to be ARy = 0.5 AR. From this, the “lift” slope is found 
similarly to Cl q above, using the same infinite-span wing lift slope ao but the new ARy, and 
multiplying by the area ratio Sy/S since this sideslip is based on the vertical tail area: 

C ° e = I + a Q / (ire ARy) ( ~S ) (28) 

(in I/rad). The tail is symmetrical for stability, so naturally there is no side force when f3 = 0, 
leading to the force equation 

C c (P) = Cc e P. (29) 
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The drag resulting from this side angle manifests like the induced drag resulting from a lift force: 


Cd, 


C ? 


c 


iveAR 


v 



(30) 


(in 1/rad 2 ). The overall aerodynamic force coefficient formulas for use in equations (19-21) are 
therefore 


Cd{(m 1 pj /f «• ) 


Cc(Pw) 
Cl ( lh w ) 


Cd 0 + C dL (C L (a w ) - C Lmin f + ( C L (a w )) 2 

+ (dw) (|) 

Cc'g Pw 

Clc («id - «o). 


( 31 ) 

(32) 

(33) 


The model developed herein is compared with that of [1] in Table 1. This table displays 
several typical gliding performance benchmarks for wings-level, trimmed, non-accelerating flight, 
both for this model and that of Beron-Rawdon’s study, for a selected set of vehicle parameter 
combinations {£, AR , W} reported in [1], The performance benchmarks are maximum lift/drag 
ratio (L/D) max , and minimum sink rate (— h) m ; n . The table also provides the speed V for these 
flight conditions with the specified parameters, and V for the case where the vehicle sink rate 
is 2(— /i) m ; n . Note that the values of (— h)min in the results from [1] remain the same for each 
{£, AR} combination. This is because, in that study, W was chosen to result in those specific 
values. As can be seen from the table, there is good general agreement between our values and 
those from [1], The principal qualitative divergence is that our values of ( L/D ) ma x do not vary 
with W, while those from [1] do, due to the latter study including Reynolds number effects which 
were ignored in this study. On the other hand, this model provides all necessary aerodynamic 
parameters for point-mass flight dynamics as continuous functions of the vehicle parameters, 
while retaining sufficient agreement with data from [1] and [8] to confidently use the model for 
exploration of small-glider flight dynamics and performance. 

In using these models, it must be stressed that the linear and quadratic angle dependencies 
in the force formulas are approximations which are accurate for small angles; in reality the lift 
will not actually increase unboundedly as a w does, for example. Thus, trajectories for this 
model which involve extreme values of a w and fi V! will be suspect, but for reasonably small 
values of the control angles these formulas should be accurate enough and appropriate to use. 
This is especially relevant in the case of smaller, slower gliders, since the aircraft velocity may 
be on the same order of magnitude as that of the arbitrary ambient wind field through which it 
flies, with the result that the wind-relative a w and values that the aircraft may transiently 
experience could be well outside the range of validity of the quadratic drag/linear lift/linear side 
force model representation. 


Section 6: Glider Model Test Simulations 


To test out the glider model equations of motion and the Fortran codes implementing them, 
we have performed simulations for some simple example problems using the codes. These were 
done to observe the behavior of the model qualitatively and see how it matches expected behavior 
for the glider. 
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Table 1: Comparison of the Current Model Flight Condition Values with Beron-Rawdon Results, 


over a Range of Glider Configurations. 
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A series of simulations were performed which calculated trajectories for gliders given con- 
stant control angle inputs at various conditions. The glider used for these simulations had values 
of i = 60 in, AR = 16, and W = 3 lb. The glider was given constant velocity-relative control 
angles of a = 10°, j3 = 10°, and a = —10°, and initial values of altitude h = 200 ft, velocity 
V = 47.9 ft/s, elevation angle 7 = —3.32°, and azimuth angle x = 0°. In the absence of ambient 
wind these conditions produce a quasi-steady-state solution, with the glider spiraling downward 
around a cylinder of constant diameter of about 800 ft, with V and 7 remaining constant. The 
trajectory is shown for 60 seconds of flight in Figure 1 (in two dimensions, from above) and 
Figure 2 (in three dimensions). 



Figure 1: Constant- Control Glider Trajectories - Top View (Configuration 1). 


With the addition of a constant-strength w y = —5 ft/s crosswind, the glider drifts off from 
that cylindrical descent, as the aircraft is aided by the wind during certain portions of the 
trajectory. This happens specifically when the glider is moving in the negative y direction, with 
the wind behind it causing it to pick up speed, which then also carries over into the negative x 
direction part of the trajectory. The glider is slowed down while flying in the positive y direction, 
impeded by the wind, leaving it slower along the trajectory moving in the positive x direction. 
In sum this results in a translational motion in the negative y direction being adding to the 
spiraling path, along with a lesser motion in the negative x direction. The glider also loses 
slightly more altitude than in the zero- wind case, as the wind makes the glider diverge from the 
more efficient quasi-steady trajectory. 

The third simulation to be performed used a wind gradient of w y = — 0.025*/r ft/s (stronger 
crosswinds at higher altitudes). At the initial altitude the wind velocity is w y = —5 ft/s, the 
same as the constant-wind case, and so those two trajectories are similar at first. Then as the 
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■ - ' - Zero Wind 

Wind w =-5 ft/s 

Wind w V =-0.025*h ft/s 

— y 


0 


Figure 2: Constant-Control Glider Trajectories (Configuration 1). 


glider loses altitude and the wind becomes slower, the trajectory’s drift slows and it approaches 
more of a steady spiral like the zero-wind case. The differences between the three trajectories 
can be seen in Figures 1 and 2 . 

The above simulations are also replicated using a different set of glider parameters: l = 100 
in, AR = 16, and W = 4.5 lb. The larger and heavier glider can catch more passing air and 
thus gain more force for maneuvering, making it easier for the wind to change the glider’s 
course. As an example of this, with zero wind the second glider configuration has a quasi-steady 
spiraling path of diameter roughly 100 ft, about one third that of the first glider configuration 
(of parameters 60/16/3). Although the maximum L/D ratio is almost the same for the two 
configuration parameter sets, the minimum sink rate is 1.56 ft/s for the second set as compared 
to 2.14 ft/s for the first. The effect of this can be seen in that the second configuration only 
loses about 80 ft of altitude during the same sixty seconds of zero- wind descent in which the first 
configuration loses about 130 ft. In this case the initial conditions for the quasi-steady spiral 
with velocity-relative control angles {a, / 3 , a} = { 10 °, 10 °, — 10 °} are altitude h = 200 ft, velocity 
V = 33.8 ft/s, elevation angle 7 = —3.04°, and azimuth angle x = 0°. 

Configuration 2 is also used with a constant wind w y = —5 ft/s, and a wind gradient 
w y = —0.025 * h ft/s (both as for configuration 1). The three trajectories are plotted in Figure 
3 (two dimensions from above) and Figure 4 (three dimensions). As with configuration 1 , the 
spiral predominantly drifts in the negative y direction, and to lesser extent in the negative x 
direction, as the wind affects the glider’s momentum. This is a result, as before, of the changing 
vehicle orientation with respect to the wind. However, it looks much more extreme here, since 
there is close to the same amount of drift but reflected on a smaller diameter spiral in the flight 
path. The wind gradient case also shows, as before, an initially similar path with drift due to 
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-500 -400 -300 -200 -100 0 


Distance y (ft) 

Figure 3: Constant-Control Glider Trajectories - Top View (Configuration 2). 



Figure 4: Constant-Control Glider Trajectories (Configuration 2). 
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the wind, before the drift slows due to the lesser winds at lower altitude. 

Figure 5 plots the zero- wind and constant -5 ft/s wind trajectories for both glider config- 
urations, to show the tighter, faster spirals of configuration 2 compared with the slower, larger 
spiral of configuration 1. As can be seen here, the glider flight behavior varies depending on the 



Figure 5: Constant-Control Glider Trajectories - Top View (Both Configurations). 


particular glider characteristics and wind strength. Even more extreme differences can result 
from larger changes in these parameters. For example, greater winds or a lighter, smaller vehicle 
can result in the glider being completely blown along with the wind, unable even to complete 
a turn back into the wind with the given control inputs. The particular control inputs used 
obviously will affect the flight behavior as well. If more extreme control inputs are used, then 
the wind will strike the aircraft from different angles and the spiral/drift patterns may be greatly 
altered. For example, a tighter spiral can result in the drift being mostly in the negative x direc- 
tion, as the accumulated speed has little chance to exert itself before the turn from negative-y 
to negative- x. A large change in the speed of descent may result as well for changes in the flight 
conditions. 

Section 7: Conclusions 

We have presented here a model of glider flight in the presence of ambient winds, which in 
simulations has resulted in realistic flight properties. The model includes a detailed accounting 
of the wind effects explicitly in its formulation, not merely as perturbations from a zero-wind 
base system of equations. Control of the system is achieved through regulation of the velocity- 
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relative aerodynamic attack, side, and bank angles, which in turn influence the strength and 
direction of the aerodynamic force vector. 

The comparison of (L/D) max , ( — /i) m ; n , and 2 ( — /i) m ; n (and the associated velocities) as 
calculated by the model described here with their values in [1] results in generally good corre- 
spondence. Even where there are differences, a qualitative behavioral match is seen, with the 
calculated values following similar patterns of change with respect to the parameters {£, AR , W } 
as the Beron-Rawdon values do. The test example sets display the flight dynamics behavior of 
the glider model in different ambient wind conditions and with different physical glider charac- 
teristics, looking at the forward solution of the flight dynamics with fixed control input. The 
generated trajectories show that the model results in believable flight behavior in three dimen- 
sions with a range of wind types and for different structural glider characteristics. 

With some relatively simple examples having been tested and the simulated glider having 
reacted in an expected manner, the model is now ready to be used in more complex and useful 
problems. Three-dimensional trajectory optimization procedures are being studied; in particular 
we are looking at “dynamic soaring” type behavior, where a maneuverable aircraft can extract 
energy from a wind gradient and thus minimize its loss of altitude (or even gain altitude) without 
thrust and in purely horizontal winds. An alternate scenario is flight through a region containing 
a number of scattered thermal updrafts, with the question of how to optimize the use of such 
thermals while at the same time following close to a desired flight path. Another avenue of 
research involves the extension of the model to consider uncertainties in the wind field or other 
elements, to more realistically represent conditions in the real world. Related to this would be 
the construction of a feedback control, based on open-loop controls generated using optimization, 
to allow the glider to respond to unexpected changes or uncertainties in the environment. 

Appendix I: Detailed Matrix Rotation Formulation 

The Euler rotation sequence {x, 7, cr} is expressed by the matrix Ryi , and moves from 
the inertial frame to the velocity frame. The x rotation about the axis z is first, followed by 
the 7 rotation about the y- axis in the new post-x frame, followed by the a rotation about the 
a:- axis in the post-x-post-7 frame. This involves the use of coordinate transformation matrices 
to ensure that the rotations are done in the correct frame of reference. Therefore, starting from 
the basic inertial axes Ui, the x- rotation matrix R x takes these to the post-x frame axes (in 
inertial coordinates): U a = R x Ui ■ The transformation matrix from these new U a coordinates to 
inertial coordinates is given by that set of axes ( Tj a = U a ). 

Since the next rotation (through 7) is done in the post-x frame, the Tj a coordinate trans- 
formation must be used before and after the rotation, resulting in U b = (T 1 aR^T^J )U a = tj a R 1 \ 
and Tib = Ub- The a- rotation similarly must have the Ti b transformation before and after it to 
change into post-x-post-7 coordinates and back. This results in Uy = (T/ b R a TJ- b l )U b = U b R a . 
Now, since overall we are in inertial frame coordinates, defined as C/j = /, this means U„ = R x 
and U b = U a R 7 = ^x^'f Therefore the final result of the sequence of three rotations in three 
different reference frames is Ryj = Uy = U b R a = R x R 7 R a . 

This shows that the combined rotation process from the inertial frame axes to the velocity 
frame axes is equivalent to doing the three individual matrix rotations all about the original 
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inertial frame axes (z, y, and x). but in the reverse order, as given by: 



cos X 

-sinx 

0 ' 


cos 7 

0 

sin 7 


' 1 

0 

0 

Rvi = 

sinx 

cos x 

0 


0 

1 

0 


0 

cos a 

— sin <7 


0 

0 

1 


— sin 7 

0 

cos 7 


0 

sin ct 

cos a 


(34) 


The rotation Rbv through —[3 and a is given by: 
Rbv = 


cos a 

0 

sin a 


cos (3 

sin /3 

0 ' 

0 

1 

0 


— sin / 3 

cos (3 

0 

— sin a 

0 

cos a 


0 

0 

1 


( 35 ) 


This is not an Euler sequence as in equation ( 34 ) where the three angles are reversed due to the 
nature of the rotation sequence. Here both rotations are considered to be about the velocity 
axes, since they represent comparatively small divergences from those axes in the horizontal and 
vertical directions, and so they are expressed appropriately in this order. 

Similarly to equations ( 34 - 35 ), the sequences of rotations through the wind- relative frame 
are given by: 



cos Xw - sin Xw 0 


cos Xw 0 sin Xw 

R\vi — 

sin Xw cos Xw 0 


0 1 0 


0 0 1 


— sin Xw 0 cos Xw 


cos a w 0 sin a w 


cos ( 3 W sin ( 3 W 0 

Rbw = 

0 1 0 


— sin / 3 W cos ( 3 W 0 


— sin a w 0 cos a w 


0 0 1 


1 

0 

0 

0 

cos a w 

— sin a w 

0 

sin a w 

cos a w 


( 36 ) 


( 37 ) 


The process of finding the unknown wind-relative angles for evaluating the equations of 
motion begins with the knowledge that x w = R\\[X. from which we can find Xw and j w . 
Componentwise this requirement is 


with x w 


x w 


COS Xw cos 7 w 
sin Xw cos j w , 
- sin x w 


V w / 1 1 V w 1 1 . The angles are computed by 
7 w = sin _1 (-(A w ) 3 ), 

Xw = sign{(x w ) 2 / cos x w )cos~ 1 {{x w ) 1 / cos x w )- 


( 38 ) 


( 39 ) 

( 40 ) 


The value of Xw found by the inverse sine is unique, due to the restriction 7,,, < tt/ 2 , which 
maintains the uniqueness of the rotational angle sequence. The sign check in equation ( 40 ) 
ensures that the correct inverse cosine solution has been used. 

To find the remaining wind- relative angles we use equation ( 18 ), which can be slightly 
rewritten (moving the known Xw and Xw terms to the opposite side of the equation) to produce: 


M 


cos Xw 0 — sin Xw 

0 1 0 
sin Xw 0 cos Xw 

1 0 0 
0 cos a w — sin a w 
0 sin a w cos a w 


cos Xw sin Xw 0 

— sin Xw cos Xw 0 

0 0 1 

cos a w 0 sin a w 
0 1 0 

— sin a w 0 cos 07, 


RviRbv 

cos (3 W sin j3 w 0 
— sin f3 w cos f3 w 0 

0 0 1 


• (41) 
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The third column of this matrix equality is given by 


^(1..3,3) 


sin a w 

— sin a w cos a w 
cos a w cos a w 


Thus we can calculate a w and a w by: 

a w = sin _i (M (li3) ), 

g w = sign(-M( 2 , 3 )/ cosa w )cos~ 1 {M (3 ^)/ cosa w ). 


(42) 


(43) 

(44) 


As with equation (39), the value of a w in equation (43) is found uniquely by the inverse sine 
function since |a w | < tt/2. Also, as in equation (40), the sign check in equation (44) serves to 
specify the correct inverse cosine solution for a w . The angle j3 w is found in a similar manner, 
starting with the first row of the matrix equality in equation (41): 


. 3 ) = [ cos cx w cos fdw 

cos a w sin / 3 W sin a w j , 

(45) 

and evaluating the equation 



Pw = sign{M (1 ^/ cos a w ] 

1 cos _1 (M( lil ) / cos a w ). 

(46) 


With the values for all of the wind-relative angles found, the aerodynamic force terms can be 
found using equations (19-23). 
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